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Abstract 

A three-dimensional unstructured grid approach to aerodynamic shape sensitivity analysis and design 
optimization has been developed and is extended to model geometrically complex configurations. The advantage of 
unstructured grids (when compared with a structured-grid approach) is their inherent ability to discretize irregularly shaped 
domains with greater efficiency and less effort. Hence, this approach is ideally suited for geometrically complex 
configurations of practical interest . In this work the nonlinear Euler equations are solved using an upwind , cell-centered, 
finite-volume scheme. The discrete , linearized systems which result from this scheme are solved iteratively by a 
preconditioned conjugate-gradient-like algorithm known as GMRES for the two-dimensional geometry and a Gauss-Seidel 
algorithm for the three-dimensional ; similar procedures are used to solve the accompanying linear aerodynamic sensitivity 
equations in incremental iterative form. As shown, this particular form of the sensitivity equation makes large-scale 
gradient-based aerodynamic optimization possible by taking advantage of memory efficient methods to construct exact 
Jacobian matrix-vector products. Simple parameterization techniques are utilized for demonstrative purposes . Once the 
surface has been deformed ’ the unstructured grid is adapted by considering the mesh as a system of interconnected springs . 
Grid sensitivities are obtained by differentiating the surface parameterization and the grid adaptation algorithms with 
ADIFOR (which is an advanced automatic- differentiation software tool). To demonstrate the ability of this procedure to 
analyze and design complex configurations of practical interest, the sensitivity analysis and shape optimization has been 
performed for a two-dimensional high-lift multielement airfoil and for a three-dimensional Boeing 747-200 aircraft. 


1. Introduction 

With the aid of modem computers, aerodynamic design 
optimization procedures [1-5] have emerged which directly 
couple the fields of computational fluid dynamics (CFD), 
sensitivity analysis, and numerical optimization. These 
procedures have enormous potential as design tools and are 
therefore receiving considerable attention in the 
aerospace, automotive, and biomedical research 
communities (among others). Bottlenecks, moreover, 
associated with the analytic evaluation of discrete 
sensitivity derivatives, appear to have been address [3] via 
the use of an incremental iterative solution of the 
sensitivity equation [5] where memory efficient methods 
[6] are used to construct Jacobian matrix-vector products. 


* Graduate Research Assistant. Student Member, AIAA. 

{ Associate Professor. 

Professor. Fellow, AIAA. 

Copyright © 1997 by James C. Newman, III. Published 
by the American Institute of Aeronautics and Astronautics, 
Inc. with permission. 


Solutions to the excessive CPU run times, to perform the 
design optimization, are being explored through the use of 
simultaneous analysis and design optimization (SAADO) 
[7], one-shot methods [8], and exploiting parallel 
computing architectures [9,10]. Another crucial hurdle, for 
these aerodynamic optimization procedures to become 
useful design tools, is their ability to analyze and design 
complex configurations of practical interest. Elliot and 
Peraire [11], with regards toward the geometrically 
complex domains associated with the integration of the 
engine into the wing design process and to the possible 
mutlipoint design of the aircraft’s high lift system and 
cruise design, assert that this may be "the step that 
determines the economic viability of the vehicle ” . 

As recently noted by Reuther et al. [2] " while flow 
analysis has matured to the extent that Navier-Stokes 
calculations are routinely carried out over very complex 
configurations, direct CFD based design is only just 
beginning to be used in the treatment of moderately 
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complex three-dimensional configurations'*. This is 
primarily due to the fact that to generate a single structured 
grid about such a configuration is difficult, if not 
impossible. Thus, to handle a typical complex geometry 
of practical interest, some sort of domain decomposition 
scheme must be incorporated into the design code. For 
structured grid solvers, these techniques would include 
multiblocked, zonally patched, and overlapped (sometimes 
referred to as Chimera) grid algorithms. However, as the 
geometric flexibility of the method increases, so does the 
complexity of the underlying algorithm. Since the use of 
sensitivity analysis, to evaluate the needed gradients for a 
numerical optimizer, is still evolving, little work has been 
done toward extending these algorithms to include these 
domain decomposition methods. The research which has 
been accomplished has mostly concentrated on the use of 
multiblocked grids. To this end, Reuther et al. [2] have 
developed a multiblock-multi grid adjoint solver 
(“continuous” or “control theory” approach [12]) which 
was applied to the wing redesign of a transonic business 
jet. Eleshaky and Baysal [13] developed a multiblock 
“discrete” adjoint solver which was applied to a simple 
axisymmetric nozzle near a flat plat. As for the use of the 
more advanced domain decomposition methods (zonal and 
overlapped grids), and combinations of the three various 
types, Taylor [14,15] has differentiated an advanced flow- 
analysis code to perform the discrete sensitivity analysis. 

An alternative, to resorting to structured grid domain 
decomposition methods to cope with complex 
configurations, is the use of unstructured grid schemes. 
Since triangles and tetrahedra are the simplest geometric 
shapes possessing area and volume, respectively, they are 
capable of resolving irregularly shaped domains easier and 
with greater efficiency. Another attribute of unstructured 
grids is that they may be adapted and locally enriched 
where needed without affecting other regions of the mesh. 

As for unstructured grid approaches to aerodynamic 
design optimization, Beux and Dervieux [16] performed 
spatially first-order accurate sensitivity analysis and 
optimization of a two-dimensional nozzle using a 
continuous adjoint method to derive the optimality 
conditions, but a discrete approach for computer 
implementation. Newman, Taylor, and Burgreen [17] 
subsequently developed a two-dimensional, and later a 
three-dimensional [3], second-order spatially accurate 
discrete sensitivity analysis approach which has been used 
to perform the design optimization of airfoils and 
transport wings in transonic flow. Elliot and Peraire [11] 
have also developed an unstructured discrete sensitivity 
analysis approach which was used to match target pressure 
distributions for a two-element airfoil, a 3D infinite wing, 
and a wing-body configuration. Subsequently, Elliot and 
Peraire [4] have applied their algorithm to perform the 
inverse pressure design of a business jet wing immersed in 
transonic flow. An equally impressive use of unstructured 
grid approaches, for the design of geometrically complex 
devices, has been performed by Burgreen and Antaki [18]. 
In Ref.[18], CFD-based design optimization methods are 


used to improve the thrombogenic performance of an axial 
flow blood pump. The research of Burgreen and Antaki 
[18], furthermore, represents the expansion of traditional 
aerodynamic design optimization procedures into the 
biomedical field to aid in artificial heart design. More 
recently, Anderson and Venkatakrishnan [19] have 
developed an unstructured grid approach to sensitivity 
analysis which truly utilizes a continuous adjoint 
approach. Moreover, in Ref.[19], limitations of the 
continuous adjoint approach are discussed and a hybrid 
continuous-discrete approach, which addresses some of 
these deficiencies, is developed. 

In this work the current unstructured grid approach to 
aerodynamic design optimization [3,17] is demonstrated 
on non-trivial, complex configurations. Presented herein 
is a discussion of the algorithms used to solve the 
nonlinear fluid and the linear sensitivity equations, to 
parameterize the design surfaces, and to perform the 
unstructured grid adaptation. Special considerations that 
arise from the use of unstructured meshes, as well as the use 
of memory efficient Jacobian matrix-vector product 
methods which make large-scale optimization possible, 
are discussed. To demonstrate this procedure, the 

aerodynamic shape sensitivity analysis and design 
optimization of a high-lift multielement airfoil and for a 
complete Boeing 747-200 aircraft is performed. Accuracy 
is accessed by comparing the analytically obtained 
sensitivity derivatives with central finite-differences. 

2. Aerodynamic Design Optimization Problem 

Aerodynamic design optimization is simply a 

constrained minimization problem which attempts to 
reduce an objective or cost function F(Q,X,fi k ) subject to 

constraints Here,g is the aerodynamic state 

vector, X is the computational mesh, and p k is the vector 

of design variables which control the shape of the 
configuration. 

A procedure to accomplish this minimization is 
obtained by linearizing the above constrained problem and 
then solving the resulting set of equations. For a gradient- 
based optimization method, such as the Method of 
Feasible Directions [20] used in the present work, frequent 
evaluations of the objective function and constraints as 
well as sensitivity gradient information are required. The 
sensitivity gradients of the objective function, VF, and 
the constraints, VCy, are commonly referred to as the 

sensitivity derivatives. These sensitivity derivatives may 
be simply evaluated by finite differences; however, this 
approach is not only computationally expensive, it has 
been found at times to produce highly inaccurate gradient 
approximations. The preferable approach is to obtain the 
sensitivity derivatives quasi -analytically via 
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To compute the sensitivity derivatives in Eqs.(la,b), the 
sensitivity of the state vector dQjdfi k is needed. This, 

consequently, results in the difficulty of solving an 
extremely large system of linear equations. The methods 
used in the present work to obtain the state vector, and the 
sensitivity of the state vector, will be discussed in 
sections to follow. 


3, Fundamental Equations 


Nonlinear Aerodynamic Analysis 
The nonlinear fluid model considered in this work will 
be the three-dimensional Euler equations. These equations 
represent the conservation of mass, momentum and energy 
for an inviscid compressible flow. Applying the backward 
Euler time-integration scheme to the unsteady term and 
linearizing the right hand side of the semi-discrete 
approximation yields 
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where /?/ represents the steady-state residual 

R i =ffEhdA = Y, E i.j (3) 
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with the face area through which the flux passes contained 
within E . 

In the present work, the inviscid flux vector, E, and the 
Jacobian, dR-JdQ, are evaluated using the flux vector 

splitting technique of Van Leer [21]. The Jacobian matrix 
may then be expressed in terms of the Van Leer fluxes as 


dR L = y !Ek jg li , d El jg LL 

dQ ji(i)[dQJ.j dQ dQjj dQ , 
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where the flux Jacobians are evaluated with variables 
interpolated to the j cell faces, and dQj jjdQ represents 

the cell center contribution from the interpolation or 
reconstruction. When higher-order spatial accuracy is 
desired in Eq.(4), the form of the Jacobian becomes 
extremely complicated and the computational stencil very 
large. This is due to the upwind biased interpolation 
scheme that must be used for unstructured grids. The order 
of accuracy of the aerodynamic analysis, however, is 
determined from the evaluation of the residual vector, and a 
first-order Jacobian has been found sufficient to converge 
the implicit scheme. Since the left hand side operator is 
not an exact linearization of the residual, the ability to 
achieve quadratic convergence is now lost. It should be 
noted that an inexact linearization is not permitted for the 
aerodynamic sensitivity equation. This is because the 
underlying equations are linear and no approximations to 


the higher-order Jacobian matrix are allowed.. This will be 
discussed in a subsequent section. 

A higher-order upwind scheme is obtained by expanding 
the cell-centered solution to the left and right of each cell 
interface using a Taylor series expansion [22], This 
expansion may be expressed as 

G/ = Q+VQ*& (5) 

where the solution gradient, Vg » at the center of the cell 

is found using the geometric invariant features of 
tetrahedra. The expression for the solution gradient at the 
cell center may be obtained from application of Green's 
theorem as 

VQ.tf = -j |j(G„| +Qn2 +G^)-G„4 (6) 

where Q nV Q n 2 , Qn$ are the primitive variables at the 

three nodes that constitute the face through which the flux 
passes, Ar is the distance from the centroid of the 
tetrahedron to the center of that face, and Q n4 are the same 
variables at the fourth node of the tetrahedron. 

The data at the nodes are obtained from the cell centered 
solution by using either an inverse distance or a psuedo- 
Laplacian weighting procedure [23]. Both procedures, 
described in Ref.[24], attempts a multidimensional 
weighted averaging of the form 



where w/ is the computed weighting factor from the desired 
node, n, to the surrounding N cell centers. 

Aerodynamic Shape Sensitivity Analysis 
As noted in a previous section, to determine the needed 
sensitivity derivatives, the sensitivity of the state vector 
dQfdfi k is required. To obtain this, the discrete residual 
vector (for a steady-state solution) may be recast as 

*{m)’X(h)’fo)= o < 8 > 

where the explicit and implicit dependencies of the residual 
on the state vector, the computational mesh, and the 
design variables are asserted. 

At this point, one of two discrete formulations may be 
used to determine the sensitivity derivatives. These 
formulations are referred to as the direct differentiation 
method and the adjoint variable method. For reasons which 
will be summarized in the conclusions, the direct approach 
is used in the current work. For a more detailed discussion 
of both methods, and their associated boundary 
conditions, the reader is referred to Ref .[25]. 

In this formulation Eq.(8) is directly differentiated with 
respect to the vector of design variables and rearranged to 
produce the following linear equation 
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dR dQ _ ( dR dR dX^ 

dQ dpk (dp k + dXdp k ) 

where dRjdQ and dRjdX are the Jacobian matrices 
evaluated at a converged flow solution, and dXjdp k is the 
grid sensitivity term. It should be noted that the task of 
constructing exactly or analytically all of the required 
Jacobians and derivatives by hand, and then building the 
software for evaluating these terms can be extremely 
complex. This problem is exemplified by the inclusion of 
even the most elementary turbulence model (for viscous 
flow) or use of a sophisticated grid generation package for 
adapting (or regenerating) the computational mesh to the 
latest design. A^ promising possible solution to this 
problem, however, has been found in the use of a technique 
known as automatic differentiation (AD). AD involves the 
application of a precompiler software tool called ADIFOR 
(Automatic Differentiation of FORtran, Ref. [26]). This 
software has been utilized, with much success, to obtain 
complex derivatives from advanced CFD and grid 
generation codes for use within aerodynamic design 
optimization procedures [14,15,27-29]. 

In the present work, the Jacobians dRjdQ y dR/dX as 
well as all derivatives (except for the grid sensitivity term) 
are constructed by hand. This is due to the fact that an 
inviscid fluid model is assumed, with the inviscid fluxes 
being constructed via the flux vector splitting scheme of 
Van Leer (a scheme which is continuously differentiable 
and well documented). ADIFOR, on the other hand, is used 
on the unstructured grid adaptation algorithm to provide 
the required grid sensitivity terms. Details of this 
algorithm and the evaluation of grid sensitivities will be 
discussed in a later section. 

The solution of Eq.(9) above poses the difficulty of 
solving an extremely large linear system for each design 
variable. Solving these systems, however, is made more 
tractable when the above equations are recast into what has 
been termed the incremental iterative form [5,14,2930] as 
follows 
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where A may now be any convenient approximation to 
the higher-order Jacobian which converges the linear 
system. This is because the equations are now cast in delta 
form, with the physics contained in the right-hand-side 
vector. It has been found that the first-order Jacobian 
works well for use in the coefficient matrix of Eq.(lOa), 
and is therefore used in the present work. A more detailed 
and thorough discussion of this incremental iterative 
technique may be found in the above cited literature. 

Two particularly attractive features of the incremental 
iterative strategy are that (0 a more diagonally dominant 


matrix may be used to drive the solution of the linear 
systems (as opposed to the sometimes ill-conditioned 
higher-order Jacobian), and (ii) the higher-order Jacobian 
now resides on the right-hand-side of the equations and 
may be dealt with in an explicit manner. When in this 
form, only the A:-vectors resulting from the matrix-vector 
product of [dRjdQ) • (dQjdfi k ) are of concern. Hence, CPU 

time and memory efficient methods for constructing the 
-exact matrix- vector product can be utilized. To this end, 
Barth and Linton [6] have developed a new technique which 
permits the construction of the higher-order Jacobian- 
vector product using slightly less memory than that which 
would be required to evaluate the first-order Jacobi an- 
vector product. This is accomplished by avoiding the need 
to assemble the full Jacobian prior to multiplication. With 
the details omitted (and the reader directed to Ref.[6] for the 
proof and further explanation of this method), this 
technique, applied to the desired matrix-vector product in 
Eq.(lOa), may be symbolically written (using the notation 
of Eq.(4)) as 


dRi dQ__ y dE U ( dQ X + 
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where ( dQ/dP k ) f j and ( dQ/dP k ) + . are the vectors 
reconstructed from dQjdP^ using the same scheme 

employed in the CFD analysis. The resulting vector, from 
this matrix-vector product, is then scattered to the adjacent 
cells in the same manner as used for the nonlinear flow- 
residual calculation. 

It should be noted that this method only requires the 
storage of the 5x5 flux Jacobians, and the reconstructed 
vectors, at the cell faces. Since this product is to be used in 
the sensitivity analysis, the memory which was utilized to 
compute the flux Jacobians for the first-order Jacobian, 
and that used to reconstruct the CFD state vector, may be 
reused. Thus, the spatially first- and higher-order accurate 
sensitivity analysis may be performed with virtually the 
same memory as the CFD analysis. The only additional 
memory is due to the storage of the grid and metric 
sensitivity terms and the derivative dRjdp k qi 
[ dRfdX)[dXjdP k ). (Note that for geometric design 
variables dR/dp k is zero, and for non-geometric design 
variables (dR/dX) [dX/ dp k ) is zero). Another attribute of 

this method is that the matrix-vector product computation 
only requires a fraction of the CPU time originally needed 
to assemble the full higher-order Jacobian; hence, the 
benefits are two-fold. The use of Barth and Linton’s 
technique within the incremental iterative method has 
significant ramifications in that it makes large-scale 
optimizations of practical three-dimensional 
configurations possible. 


4. Solution Methodology 

The linear systems resulting from the nonlinear 
aerodynamic analysis, as well as those from the 
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aerodynamic shape sensitivity analysis, can be expressed 
in the simple form Ax=b. Within the optimization 
process, it is evident that the aerodynamic analysis not 
only consumes more CPU time (than the shape sensitivity 
analysis) to converge the nonlinear systems, it also is 
needed more frequently. Thus, solution algorithms which 
have high convergence rates are imperative when multiple 
analyses are to be performed. To this end, it has been found 
that a fully implicit iterative solver, utilizing 
preconditioned GMRES techniques developed by Saad and 
Schultz [31], often out perform conventional explicit as 
well as classical iterative solvers [32-34]. In previous 
work [3,17], all linear systems resulting from the solution 
of Eq.(2) and Eq.(lOa) for the aerodynamic analysis and 
shape sensitivity analysis, respectively, utilized GMRES. 
However, the selection of the preconditioner used in 
conjunction with GMRES essentially governs the 
performance and memory required by this algorithm. 
Computations performed in Refs. [3, 17] used an 
incomplete LJU factorization (ILU(O)) as a left 
preconditioner for Eq.(2) and as a right preconditioner for 
Eq.(lOa). In the current work, the two-dimensional 
computations about the multielement airfoil use GMRES, 
but due to the memory requirements associated with this 
algorithm, it was not possible to utilize it for the Boeing 
747-200 configuration. Hence, a Gauss-Seidel iterative 
method was used to solve the fluid and sensitivity 
equations for this geometry. 

As a final note, it has been observed that the ordering of 
the cells in a grid has an affect on the rate of convergence 
of iterative solvers [35]. With this in mind, many 
researchers [343637] currently working with unstructured 
grids (which usually have a random ordering) have adopted 
renumbering algorithms such as Cuthill-McKee (CM) [38] 
or reverse-CM [39]. These algorithms attempt to reorder a 
given mesh such that the bandwidth of the coefficient 
matrix is minimized. In the present work, reordering is 
accomplished, during the initial preprocessing of the grid, 
using the Gibbs- Pool e-Stockmeyer [40] algorithm. 

5. Design Surface and Grid Representation 

A key aspect in any design optimization procedure is 
how the design surface and computational mesh are to be 
represented. This selection will ultimately determine (*) 
the type and number of design variables used, ( ii ) the grid 
adaptation or regeneration method, and (iii) the means 
through which grid sensitivities are calculated. In the 
following sections, the techniques used in the present 
work for (i), (ii), and (iii) will be discussed. 

Design Surface Parameterization 

Once an aerodynamic shape optimization code has been 
developed and verified, only the design surface 
parameterization routines change from application to 
application. The particular parameterization technique 
utilized depends on the geometry being studied and the 
design problem formulation. An excellent example of a 
sophisticated wing parameterization method capable of 



(a) Unstructured mesh (7614 nodes and 14919 cells). 


Upper Surface Bezier Control Point 
I Design Variable Number 3 



Design Variable Number 8 


(b) Design surface parameterization for the vane. 



(c) Grid sensitivity dy / d/i 7 . 


Fig. 1: Mesh, vane parameterization, and grid sensitivity 
for the multielement airfoil. 

modeling wing-section (airfoil) definitions, taper 
distribution, sweep, span and spanwise bending, global 
angle-of-attack, and twist schedule was developed in 
Ref.[41], and discussed at length in Ref.[l], In the current 
work, however, simpler parameterization techniques have 
been used for demonstrative purposes and to keep the 
number of design variables to a minimum for the Boeing 
747-200 configuration. 
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(a) Unstructured mesh (63828 nodes and 352547 cells). 



(b) Grid sensitivity dz/dfij. 

Fig. 2: Surface mesh and grid sensitivity for the Boeing 
747-200 configuration. 

For the high-lift multielement airfoil, the upper and 
lower surfaces of the vane were parameterized with separate 
Bezier curves. Details of this type of parameterization for 
the design of airfoils may be found in Refs. [ 17,32]. The 
design variables are the vertical locations of the 10 
interior Bezier control points. The unstructured mesh, and 
this parameterization, are shown in Fig. la and lb, 
respectively. As for the Boeing 747-200 geometry, the 
dihedral and twist schedule along the wing, outboard of the 
outermost engine nacelle, was parameterized with cubic 
polynomials. At the point where the engine strut meets the 
wing, point and slope continuity are enforced. Thus, only 
the coefficients of the quadratic and cubic terms are free and 
are therefore chosen as the design variables. This 
constitutes 2 design variables for the dihedral and 2 for the 
twist. The unstructured surface mesh for the Boeing 747- 
200 configuration is depicted in Fig. 2a, and is derived 
from the model tested in the NASA Ames 1 1 foot Transonic 
Pressure Tunnel (Test AR0502). It should be noted that the 
twist distribution of this geometry is not representative of 
the normal shape of the production wing. 


Unstructured Grid Adaptation 

The mesh movement strategy adopted considers the 
mesh as a system of interconnecting springs. This system 
is constructed by representing each edge of each 
tetrahedron by a tension spring. In the current method, the 
stiffness of the springs is assumed inversely proportional 
to the length of its edge. Then, for each mesh point, the 
external forces due to the connecting boundary springs are 
summed and resolved into Cartesian components. The 
result is a set of linear systems which may be solved for 
the displacements of each node using several Jacobi 
iterations. For further details on this method the reader is 
directed to the literature [42-44], Reference [44] has the 
added advantage of edge (2D; face in 3D) swapping based 
on the Delaunay criterion which greatly improves the 
performance of the method. This technique, however, is 
not currently implemented in the present grid adaptation 
algorithm, but is deemed a vital improvement. 

Grid Sensitivities 

Efficient and accurate evaluation of grid sensitivities is 
an extremely important and vital aspect in any design 
optimization procedure (which uses discrete sensitivity 
analysis). The technique used to obtain the grid 
sensitivities from the unstructured grid adaptation 
procedure results from the direct application of ADIFOR. 
Here, the subroutines which define the design variables 
P k , and the subroutines which perform the unstructured 
grid adaptation to produce the mesh X(Pk)* are 

differentiated using ADIFOR. The result is an additional set 
of subroutines which, upon compilation and execution, 
will return the grid sensitivities, dX/dp k . 

To verify that these sensitivities were indeed correct for 
both geometries, the design variables were perturbed, the 
grid adapted, and the grid sensitivities calculated via 
ADIFOR generated subroutines. These sensitivities were 
then compared with those obtained using finite -difference. 
Quantitatively, ADIFOR generated grid sensitivities 
matched finite-difference to approximately 8 significant 
digits. A qualitative representation of computed grid 
sensitivities are depicted in Fig. lc and Fig. 2b for the 
multielement airfoil and Boeing 747-200 configurations, 
respectively. Figure lc illustrates the sensitivity of the 
internal mesh to one of the design variables (Bezier 
control points) on the lower surface of the vane. Figure 2b 
represents the sensitivity of the surface grid to a twist 
design variable (specifically the coefficient of the 
quadratic term in the spanwise twist schedule). 

6. Results and Discussion 

In the present work, aerodynamic shape sensitivity 
analysis and design optimization for two sample complex 
configurations are examined: a two-dimensional high-lift 
multielement airfoil and a three-dimensional Boeing 747- 
200 aircraft. The flow conditions for the multielement 
airfoil calculation were a subsonic Mach number of 0.20 at 
16.02 degrees angle-of-attack. For the Boeing 747- 


6 




Fig. 3: Initial and optimized vane shapes. 



Fig. 4: Pressure coefficient distributions on the 
multielement airfoil. 

200, a transonic Mach number of 0.84 was chosen with a 
freestream angle-of-attack of 2.73 degrees. To verify the 
accuracy of present discrete sensitivity analysis approach, 
the sensitivity derivatives of lift-to-drag ratio, with 
respect to the geometric design variables previously 
discussed, are computed and compared with finite- 
differences. It should be noted that the work associated 
with computing sensitivity derivatives via the direct 
differentiation method does not scale with the number of 
output functions. Hence, for the computations shown, any 
number of output function sensitivities (i.e., for lift, drag, 
pitching moments.. .etc.) can be computed with little 
effort. To evaluate the corresponding finite-difference 
derivatives, however, requires two analysis runs per design 
variable for central differences. Thus, due to the expense of 
finite-difference derivatives, only those derivative 
comparisons with respect to one of the geometric design 
variables were performed. 

For the multielement airfoil, the geometric design 
variable selected was the same Bezier control point used 


above in Fig. lc for the grid sensitivity illustration. For 
this configuration the analytically obtained derivative 
dC L /dfi 7 was computed to be 0.210754 as compared to 
the central finite-difference value, 0.210733. The 
geometric design variable selected for the Boeing 747-200 
was, once again, the one presented above in the grid 
sensitivity demonstration of Fig. 2b; namely, the 
coefficient of the quadratic term in the spanwise twist 
schedule. Here, the computed derivative d(L/D)/dfii 
yielded 1.3231, and the finite-difference calculation 
1.3229. As should be expected with consistent, discrete 
sensitivity analysis, computed derivatives match finite- 
difference to approximately 4 significant digits for both 
geometries. 

Before presenting the design results, some brief words 
about the importance of design problem formulation need 
to be asserted. Sensitivity analysis is merely an extra level 
of computation that provides additional information to the 
designer. When the sensitivity analysis routines are 
coupled with the fluid solver, a mesh movement strategy, 
and a numerical optimizer, a functional design tool is 
produced. The eventual designs created with this tool will 
be only as good as the formulated design problem. If 
improperly formulated, designs can be produced that will 
violate constraints such as those needed for 

manufacturability or structural feasibility, or produce a 
design that has superb performance at one operating 
condition, but is unacceptable off the design point. Thus, 
experience is required in formulating meaningful design 
problems. 

For the multielement airfoil case, it is recognized that 
the goal of a flap system is to maintain the highest 
possible lift-to-drag ratio at the maximum lift coefficient 
[45]. With this in mind, the problem formulation 
consisted of maximizing the lift coefficient. Note, it is not 
asserted that inviscid flow analysis and sensitivity 
analysis are capable of modeling the physics or properly 
designing such a configuration (especially at high angle- 
of-attack situations such as take off and landing), but 
rather that an unstructured grid method easily discretizes 
the domain and that it is possible to carry out this type of 
design study. Furthermore, a geometric constraint has been 
placed on the thickness of the trailing edge to keep it from 
becoming excessively thick or thin. 

Results of this design study are summarized in Table 1 . 
This optimization required 7 design cycles with 83 CFD 
analyses along the line searches and a total run time of a 
little over 40.5 min. on a CRAY Y/MP. The objective 
function was increased by 6.2%. Figure 3 depicts the 
initial and final optimized vane produced by this shape 
optimization procedure. Illustrated in Fig. 4 are the 
corresponding pressure coefficient distributions about the 
multielement airfoil. It should be noted that the horizontal 
location of the vane and flap have been altered so that their 
C P distributions may be easily distinguished. As seen, the 
pressure distributions about the leading edge slat and the 
main airfoil remain roughly unchanged, but those on the 
vane and flap have been greatly altered. Another 
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(a) Initial and optimized twist and dihedral distributions. 



(b) Upstream view of the surface mesh on the outboard 
wing stations. 


Fig. 5: Optimization results for the Boeing 747-200 wing 
redesign. 

interesting design problem, which may have produced 
more dramatic improvements, could have been the design 
optimization of the shape, location, and orientation of the 
leading edge slat. 

For the Boeing 747-200 aircraft, the design problem 
was formulated to maximize the Iift-to-drag ratio. Once 
again, constraints that the lift coefficient at the final 
design be greater than the initial value, and that the drag 
coefficient be reduced for the optimized shape, have been 
incorporated. The results of this optimization are also 
shown in Table 1. It should be noted that since this design 
study was for purely demonstrative purposes, the 
optimization was halted after 3 design cycles and not 
allowed to continue until an isolated local minimum was 
found. Observe that the objective function has been 
improved by 2.7%, but at the cost of 23.4 CRAY Y/MP 
hours. This represents about 6 converged CFD analyses. 
As noted in the introduction, however, techniques are now 
being studied to reduce these excessive CPU times. 
Nevertheless, the initial and final twist and dihedral 
schedules are shown in Fig. 5a. Note that a positive twist 
angle is defined herein as leading edge up. As seen, this 
optimized wing has a greater twist at the tip station and an 


altered dihedral distribution. The initial and optimized 
surface meshes are viewed from an upstream vantage point 
in Fig 5b to show the dihedral distributions of each wing. 
Surface pressure contours for both the initial and final 
designs are illustrated in Fig. 6 and 7 for the upper and 
lower surfaces, respectively. It can be observed from close 
inspection that the upper surface of the optimized wing has 
a greater region of lower pressure than the initial wing and 
that the lower surface has a slightly higher pressure. Once 
again, a design problem which possibly could have 
produced more significant increases in the lift-to-drag ratio 
would have been to perform the shape optimization of the 
wing airfoil sections. However, this current procedure has 
demonstrated the ease with which an unstructured grid 
approach to aerodynamic shape sensitivity analysis and 
optimization may be used to analyze and design 
geometrically complex configurations. 

7. Conclusions 

A three-dimensional unstructured grid approach to 
aerodynamic shape sensitivity analysis and optimization 
has been demonstrated on geometrically complex 
configurations of practical interest. It was shown that 
shortcomings of discrete sensitivity analysis can be 
bypassed by the use of an efficient Jacobian matrix-vector 
product technique within the incremental iterative form of 
the sensitivity equation and through the use of automatic 
differentiation to obtain grid sensitivities. This approach 
was demonstrated through the shape optimization of the 
vane in a two-dimensional multielement airfoil 
configuration and by the twist and dihedral schedule on the 
outboard stations of the wing on a Boeing 747-200 
aircraft. The complexity of the geometries studied herein 
illustrate the advantages of using unstructured grids for 
aerodynamic shape optimization. 

In the current work the direct differentiation approach, 
as opposed to the adjoint variable approach, was used to 
perform the discrete sensitivity analysis. As is well 
known, for design problem formulations in which the 
number of design variables exceed the number of 
constraints plus one (for the objective function), the 
adjoint method is the preferred approach. However, when 
the reciprocal is true and there are more constraints than 
design variables, as in the case of multidisciplinary 
optimization, the direct approach is more attractive. Since 
the ultimate goal of the present work is the development of 
a multidisciplinary analysis and optimization procedure 
(using discrete sensitivity analysis), the direct 
differentiation method was adopted. Moreover, the current 
algorithm has been extended to incorporate an existing 
finite-element code to perform the aeroelastic analysis 
[46,471, and work is currently underway to perform the 
discrete aeroelastic sensitivity analysis. Results are 
forthcoming. 
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Table 1: Summary of optimization results. 



Initial 

Objective 

Final 

Objective 

Function 

Evaluations 

Gradient 

Evaluations 

Memory 1 

[MW] 

WmS^M 

MEA 

3.792 

4.028 

82 

7 

2.67/4.18 

0.062/0.675 

B747 

14.136 

14.522 

24 

3‘ 

63.3/104.1 

3.7/23.4 


* Stopped after third design cycle. 

t Memory for CFD analysis/ memory for sensitivity analysis. 

$ CPU time for converged CFD analysis/total optimization run time. 




(b) Pressure contours on the optimized geometry. (b) Pressure contours on the optimized geometry. 


Fig. 6: Upper surface pressure contours for the initial and 
optimized B747-200 configuration. 


Fig. 7: Lower surface pressure contours for the initial and 
optimized B747-200 configuration. 
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